library(tidyverse)
library(readxl)
library(GGally)
library(VIM)
library(FactoMineR)
library(stats)
library(factoextra)
library(ggrepel)Evaluating Short-Term Restoration Approaches for Degraded Grasslands in Kenya
Summary
This project evaluates the short-term effectiveness of restoration methods for degraded grasslands in Kenya. The treatments include Tilage (T), the addition of native soil inoculum (S) and Manure (M). “A grazing exclosure was installed at the campus of the Kabianga University in Kenya in April 2022, when restoration treatments were laid out and baseline soil samples taken and analysed for key soil physical, chemical and biological properties. In November 2022, fresh soil samples were taken to track short- term changes in the soil. In addition, we assessed the plant community species composition, and measured in situ soil respiration rates.” Refer “Grassland restoration experiment Kabianga Kenya.pdf” for more info.
Abrevations used
| Variables and abbreviations | Explanation |
|---|---|
| M | Manure |
| S | Soil inoculum |
| T | Tillage |
| C | Control |
| YIELD | Yield in dry weight (kg/625 m^2) |
| DM | Dry matter (%) |
| Tsoil | Temperature of the topsoil (0C) |
| Tair | Temperature of the air (0C) |
| MOIS | Moisture of the topsoil (%) |
| pH | pH in water (dimensionless) |
| EC | Electrical Conductivity (µs m-1) |
| MBC | Microbial Biomass Carbon (mg C kg soil-1) |
| MBN | Microbial Biomass Nitrogen (mg N kg soil-1) |
| Nconc | Soil total Nitrogen (%) |
| Cconc | Soil total Carbon (%) |
| P | Olsen P (μg P g-1 soil) |
| BD | Bulk density (g/cm3) |
| NH4 | Ammonium-Nitrogen (µg NH4-N g-1 DW) |
| NO3 | Nitrate-Nitrogen (µg NO3-N g-1 DW) |
Loading required packages
Loading all datasets
Some unwanted columns have also been removed and -1 is replaced as 0 for easy interpretability.
#Kabianga
kab.nov.22 <- read_excel("data/DSAS_REDEAL_Assign_Data_Aug2.xlsx",
sheet = 'Kab_nov_2022',
na = c("NA", NA)) |>
select(-c(...21,...22)) |>
mutate(M = if_else(M == -1, 0, 1),
S = if_else(S == -1, 0, 1),
T = if_else(T == -1, 0, 1))
kab.dec.23 <- read_excel("data/DSAS_REDEAL_Assign_Data_Aug2.xlsx",
sheet = 'Kab_dec_2023',
na = c("NA", NA)) |>
mutate(M = if_else(M == -1, 0, 1),
S = if_else(S == -1, 0, 1),
T = if_else(T == -1, 0, 1))
#Thurgem
thur.jun.22 <- read_excel("data/DSAS_REDEAL_Assign_Data_Aug2.xlsx" ,
sheet = 'Thur_june_2022',
na = c("NA", NA)) |>
select(-c(...13,...14)) |>
mutate(M = if_else(M == -1, 0, 1),
S = if_else(S == -1, 0, 1),
T = if_else(T == -1, 0, 1))
thur.july.23 <- read_excel("data/DSAS_REDEAL_Assign_Data_Aug2.xlsx",
sheet ='Thur_july_2023',
na = c("NA", NA)) |>
select(-c(...15,...16)) |>
mutate(M = if_else(M == -1, 0, 1),
S = if_else(S == -1, 0, 1),
T = if_else(T == -1, 0, 1))
thur.dec.23 <- read_excel("data/DSAS_REDEAL_Assign_Data_Aug2.xlsx" ,
sheet = 'Thur_dec_2023',
na = c("NA", NA)) |>
mutate(M = if_else(M == -1, 0, 1),
S = if_else(S == -1, 0, 1),
T = if_else(T == -1, 0, 1))
#Kapsarok
kaps.jun.22 <- read_excel("data/DSAS_REDEAL_Assign_Data_Aug2.xlsx",
sheet = 'Kaps_june_2022',
na = c("NA", NA)) |>
select(-c(...13,...14)) |>
mutate(M = if_else(M == -1, 0, 1),
S = if_else(S == -1, 0, 1),
T = if_else(T == -1, 0, 1))Data exploration
Missing value check
aggr(kab.nov.22, prop = F, numbers = T) aggr(kab.dec.23, prop = F, numbers = T)aggr(thur.jun.22, prop = F, numbers = T)aggr(thur.july.23, prop = F, numbers = T)aggr(thur.dec.23, prop = F, numbers = T)aggr(kaps.jun.22, prop = F, numbers = T)For the dataset kab.nov.23 the plot indicate the last row has all NA entries which is removed. And for the dataset kaps.jun.22 there are 2 NA entries in the YIELD column, here the average based on the type of treatment for which the value is missing is taken and NA was replaced. Since there are only 2 NA values, code is not written for the entire column.
kab.nov.22 <- kab.nov.22|>
na.omit()
mean_yield_1 <- kaps.jun.22 |>
filter(Code == "S") |>
pull(YIELD) |>
mean(na.rm = TRUE)
mean_yield_2 <- kaps.jun.22 |>
filter(Code == "ST") |>
pull(YIELD) |>
mean(na.rm = TRUE)
kaps.jun.22 <- kaps.jun.22 |>
mutate(YIELD = if_else(Code == "S" & is.na(YIELD), mean_yield_1, YIELD),
YIELD = if_else(Code == "ST" & is.na(YIELD), mean_yield_2, YIELD)) Checking data types
For the variables Plot No., YIELD, DM, SPE, Tsoil, Tair, MOIS, RES, pH, EC, MBC, MBN, Nconc, Cconc, P, BD, NH4, NO3, Block and for each treatment type M, S and T numeric values are expected as entry.
glimpse(kab.nov.22)Rows: 32
Columns: 20
$ `Plot No.` <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, …
$ Block <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3,…
$ Treatment <chr> "Soil", "Manure+Soil+Tillage", "Control", "Manure", "Manure…
$ Code <chr> "S", "MST", "C", "M", "MT", "MS", "ST", "T", "MST", "ST", "…
$ M <dbl> 0, 1, 0, 1, 1, 1, 0, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 1,…
$ S <dbl> 1, 1, 0, 0, 0, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0,…
$ T <dbl> 0, 1, 0, 0, 1, 0, 1, 1, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 1,…
$ Tsoil <dbl> 18.9, 18.6, 18.7, 18.7, 18.7, 18.7, 18.8, 19.7, 20.7, 21.8,…
$ Tair <dbl> 19.7, 19.8, 20.2, 20.3, 20.6, 20.9, 21.7, 24.6, 29.9, 31.2,…
$ MOIS <dbl> 53.6, 37.0, 51.7, 51.6, 40.4, 52.4, 37.6, 36.6, 40.4, 38.7,…
$ YIELD <dbl> 21.2, 13.8, 15.2, 17.8, 16.0, 16.4, 18.1, 17.1, 16.8, 15.2,…
$ NH4 <dbl> 0.7, 1.1, 2.0, 1.6, 1.2, 1.8, 0.5, 0.7, 2.1, 0.5, 1.1, 1.2,…
$ NO3 <dbl> 22.8, 31.0, 23.6, 5.3, 12.6, 24.6, 25.1, 25.3, 45.6, 14.6, …
$ P <dbl> 1.4, 2.5, 0.6, 1.4, 1.7, 1.3, 1.0, 2.1, 1.4, 1.7, 2.0, 3.2,…
$ pH <dbl> 5.0, 4.9, 4.7, 4.5, 4.6, 4.7, 4.7, 4.7, 4.6, 4.6, 4.4, 4.6,…
$ EC <dbl> 41.9, 47.4, 35.1, 11.6, 19.7, 38.5, 35.2, 36.8, 59.9, 22.0,…
$ MBC <dbl> 461.6, 512.7, 577.2, 574.4, 440.7, 514.2, 440.1, 422.5, 514…
$ MBN <dbl> 88.0, 109.9, 108.7, 109.4, 85.9, 93.4, 81.8, 83.8, 98.5, 93…
$ Nconc <dbl> 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.4, 0.5, 0.5, 0.5, 0.4,…
$ Cconc <dbl> 5.1, 5.4, 5.4, 5.0, 5.1, 5.3, 4.9, 4.8, 5.2, 5.3, 5.8, 4.6,…
glimpse(kab.dec.23)Rows: 32
Columns: 10
$ `Plot No.` <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, …
$ Block <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3,…
$ Treatment <chr> "Soil", "Manure+Soil+Tillage", "Control", "Manure", "Manure…
$ Code <chr> "S", "MST", "C", "M", "MT", "MS", "ST", "T", "MST", "ST", "…
$ M <dbl> 0, 1, 0, 1, 1, 1, 0, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 1,…
$ S <dbl> 1, 1, 0, 0, 0, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0,…
$ T <dbl> 0, 1, 0, 0, 1, 0, 1, 1, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 1,…
$ BD <dbl> 0.3589024, 0.4051220, 0.3829878, 0.4227439, 0.4143293, 0.38…
$ P <dbl> 5.9156125, 3.1895944, 1.5363629, 1.3363361, 1.8315103, 1.40…
$ YIELD <dbl> 10.43400, 6.74075, 6.88250, 6.60400, 5.60350, 6.98075, 11.3…
glimpse(thur.jun.22)Rows: 32
Columns: 12
$ `Plot No.` <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, …
$ Block <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3,…
$ Treatment <chr> "Soil", "Tillage + Soil + Manure", "Control", "Manure", "Ti…
$ Code <chr> "S", "MST", "C", "M", "MT", "MS", "ST", "T", "MST", "ST", "…
$ M <dbl> 0, 1, 0, 1, 1, 1, 0, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 1,…
$ S <dbl> 1, 1, 0, 0, 0, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0,…
$ T <dbl> 0, 1, 0, 0, 1, 0, 1, 1, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 1,…
$ Tsoil <dbl> 21.8, 22.5, 22.8, 23.5, 24.0, 24.4, 24.8, 25.7, 26.6, 27.7,…
$ Tair <dbl> 24.1, 25.7, 26.7, 27.4, 28.0, 27.8, 28.7, 30.7, 32.8, 33.8,…
$ MOIS <dbl> 44.7, 46.3, 45.0, 42.7, 28.0, 43.2, 33.3, 38.1, 33.0, 39.5,…
$ YIELD <dbl> 4.2, 5.8, 3.8, 3.1, 3.5, 5.1, 3.9, 2.5, 3.9, 4.2, 1.6, 6.7,…
$ EC <dbl> 136.0, 149.5, 165.3, 182.6, 176.0, 190.4, 152.0, 148.0, 128…
glimpse(thur.july.23)Rows: 32
Columns: 14
$ Plot_No <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 1…
$ Block <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, …
$ Treatment <chr> "Soil", "Tillage + Soil + Manure", "Control", "Manure", "Til…
$ Code <chr> "S", "MST", "C", "M", "MT", "MS", "ST", "T", "MST", "ST", "M…
$ M <dbl> 0, 1, 0, 1, 1, 1, 0, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 1, …
$ S <dbl> 1, 1, 0, 0, 0, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, …
$ T <dbl> 0, 1, 0, 0, 1, 0, 1, 1, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 1, …
$ Tsoil <dbl> 23.7, 23.6, 23.9, 24.8, 25.1, 25.0, 25.2, 26.7, 26.4, 28.0, …
$ Tair <dbl> 25.3, 26.3, 28.0, 30.4, 32.3, 32.1, 32.4, 35.3, 33.5, 33.9, …
$ MOIS <dbl> 48.1, 52.3, 47.3, 40.1, 47.9, 39.6, 47.8, 45.9, 51.6, 44.4, …
$ DW <dbl> 162.7, 214.9, 270.3, 262.5, 298.1, 648.9, 646.1, 465.1, 253.…
$ YIELD <dbl> 4.1, 5.4, 6.8, 6.6, 7.5, 16.2, 16.2, 11.6, 6.3, 8.4, 4.0, 19…
$ Nconc <dbl> 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, …
$ Cconc <dbl> 1.0, 1.1, 1.3, 1.0, 1.0, 1.4, 1.2, 1.2, 0.9, 0.9, 0.8, 1.4, …
glimpse(thur.dec.23)Rows: 32
Columns: 10
$ `Plot No.` <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, …
$ Block <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3,…
$ Treatment <chr> "Soil", "Tillage + Soil + Manure", "Control", "Manure", "Ti…
$ Code <chr> "S", "MST", "C", "M", "MT", "MS", "ST", "T", "MST", "ST", "…
$ M <dbl> 0, 1, 0, 1, 1, 1, 0, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 1,…
$ S <dbl> 1, 1, 0, 0, 0, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0,…
$ T <dbl> 0, 1, 0, 0, 1, 0, 1, 1, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 1,…
$ MOIS <dbl> 55.30603, 63.17742, 60.24427, 56.02128, 59.56154, 60.94444,…
$ P <dbl> 2.7800613, 1.6521765, 1.6484200, 1.5918039, 1.6554478, 1.41…
$ YIELD <dbl> 2.59225, 2.85375, 1.30200, 1.55025, 2.62850, 4.39375, 4.625…
glimpse(kaps.jun.22)Rows: 32
Columns: 12
$ Plot_No. <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 1…
$ Block <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, …
$ Treatment <chr> "Soil", "Tillage + Soil + Manure", "Control", "Manure", "Til…
$ Code <chr> "S", "MST", "C", "M", "MT", "MS", "ST", "T", "MST", "ST", "M…
$ M <dbl> 0, 1, 0, 1, 1, 1, 0, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 1, …
$ S <dbl> 1, 1, 0, 0, 0, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, …
$ T <dbl> 0, 1, 0, 0, 1, 0, 1, 1, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 1, …
$ Tsoil <dbl> 21.9, 21.9, 22.4, 23.5, 24.2, 24.7, 25.4, 26.2, 26.7, 28.7, …
$ Tair <dbl> 24.3, 24.9, 26.1, 27.0, 28.8, 30.9, 32.5, 31.4, 32.0, 35.6, …
$ MOIS <dbl> 12.0, 15.1, 19.9, 19.3, 25.4, 28.8, 26.2, 22.4, 18.9, 24.3, …
$ YIELD <dbl> 1.866667, 2.900000, 1.400000, 2.300000, 10.300000, 6.500000,…
$ EC <dbl> 106.0, 136.6, 140.8, 156.4, 110.8, 127.8, 101.8, 184.4, 140.…
All results were as expected.
Check for multicolinearity
Kabianga November 2022
kab.nov.22 |> select(-`Plot No.`, - Treatment, -M, -S, -T) |> ggpairs()For the dataset a high correlation is observed between,
Var 1 Var 2 r Tair Tsoil 0.909 Nconc Tair 0.615 YIELD MOIS 0.571 MBN NH4 0.651 EC NO3 0.980 Nconc NO3 0.603 MBN MBC 0.891 Cconc MBC 0.627 Cconc Nconc 0.731
Kabianga December 2023
kab.dec.23 |> select(-`Plot No.`, - Treatment, -M, -S, -T) |> ggpairs()No significant correlations were seen for the dataset.
Thurgem June 2022
thur.jun.22 |> select(-`Plot No.`, - Treatment, -M, -S, -T) |> ggpairs()For the dataset a high correlation is observed between,
Var 1 Var 2 r Tair Tsoil 0.974
Thurgem July 2023
thur.july.23 |> select(-Plot_No, - Treatment, -M, -S, -T) |> ggpairs()For the dataset a high correlation is observed between,
Var 1 Var 2 r Tair Tsoil 0.912 YIELD DW 1
Thurgem Dec 2023
thur.dec.23 |> select(- `Plot No.`, - Treatment, -M, -S, -T) |> ggpairs()No significant correlations were observed between variables.
Kapsorok June 2022
kaps.jun.22 |> select(-Plot_No., - Treatment, -M, -S, -T) |> ggpairs()For the dataset a high correlation is observed between,
Var 1 Var 2 r Tair Tsoil 0.962
Correlation between Tair and Tsoil is consistent across all datasets.
Histogram of YIELD by Treatment faceted by Block
Kabianga November 2022
ggplot(kab.nov.22, aes(x = Code, y = YIELD)) + geom_bar(color = "black", alpha = 0.7, stat = "identity") + labs(title = "Histogram of YIELD by Treaetment faceted by Block (Kabianga Nov 2022)", x = "Treatment", y = "Yield") + facet_wrap(~Block) + theme_minimal()
Kabianga December 2023
ggplot(kab.dec.23, aes(x = Code, y = YIELD)) + geom_bar(color = "black", alpha = 0.7, stat = "identity") + labs(title = "Histogram of YIELD by Treaetment faceted by Block (Kabianga Dec 2023)", x = "Treatment", y = "Yield") + facet_wrap(~Block) + theme_minimal()
Thurgem June 2022
ggplot(thur.jun.22, aes(x = Code, y = YIELD)) + geom_bar(color = "black", alpha = 0.7, stat = "identity") + labs(title = "Histogram of YIELD by Treaetment faceted by Block (Thurgem June 2022)", x = "Treatment", y = "Yield") + facet_wrap(~Block) + theme_minimal()
Thurgem July 2023
ggplot(thur.july.23, aes(x = Code, y = YIELD)) + geom_bar(color = "black", alpha = 0.7, stat = "identity") + labs(title = "Histogram of YIELD by Treaetment faceted by Block (Thurgem July 2023)", x = "Treatment", y = "Yield") + facet_wrap(~Block) + theme_minimal()
Thurgem Dec 2023
ggplot(thur.dec.23, aes(x = Code, y = YIELD)) + geom_bar(color = "black", alpha = 0.7, stat = "identity") + labs(title = "Histogram of YIELD by Treaetment faceted by Block (Thurgem Dec 2023)", x = "Treatment", y = "Yield") + facet_wrap(~Block) + theme_minimal()
Kapsorok June 2022
ggplot(kaps.jun.22, aes(x = Code, y = YIELD)) + geom_bar(color = "black", alpha = 0.7, stat = "identity") + labs(title = "Bar chart of YIELD by Treaetment faceted by Block (Kapsorok June 2022)", x = "Treatment", y = "Yield") + facet_wrap(~Block) + theme_minimal()
For all dataset the variation of YIELD across treatment across each block seem to vary slightly however it’s statistical significance need to be checked for each dataset.
Average yield across treatment
Kabianga November 2022
kab.nov.22 |> summarise(m_y_t = mean(YIELD), .by = Code) |> ggplot(aes(x = Code, y = m_y_t)) + geom_bar(color = "black", alpha = 0.7, stat = "identity") + labs(title = "Bar chart of mean YIELD by Treaetment (Kabianga November 2022)", x = "Treatment", y = "Mean Yield") + theme_minimal()
For the dataset, it can be observed that for treatments MST, MT, ST and T there is a reduction in YIELD compared to C (Control).
Kabianga December 2023
kab.dec.23 |> summarise(m_y_t = mean(YIELD), .by = Code) |> ggplot(aes(x = Code, y = m_y_t)) + geom_bar(color = "black", alpha = 0.7, stat = "identity") + labs(title = "Bar chart of mean YIELD by Treaetment (Kabianga December 2023)", x = "Treatment", y = "Mean Yield") + theme_minimal()
For the dataset, it can be observed that for treatments MSTand M there is a reduction in YIELD compared to C (Control) and T the increase is fairly low.
Thurgem June 2022
thur.jun.22 |> summarise(m_y_t = mean(YIELD), .by = Code) |> ggplot(aes(x = Code, y = m_y_t)) + geom_bar(color = "black", alpha = 0.7, stat = "identity") + labs(title = "Bar chart of mean YIELD by Treaetment (Thurgem June 2022)", x = "Treatment", y = "Mean Yield") + theme_minimal()
For the dataset, it can be observed that for treatments M, MS, MST, S, MT, ST and T there is a reduction in YIELD compared to C (Control).
Thurgem July 2023
thur.july.23 |> summarise(m_y_t = mean(YIELD), .by = Code) |> ggplot(aes(x = Code, y = m_y_t)) + geom_bar(color = "black", alpha = 0.7, stat = "identity") + labs(title = "Bar chart of mean YIELD by Treaetment (Thurgem July 2023)", x = "Treatment", y = "Mean Yield") + theme_minimal()
For the dataset, it can be observed that for treatments M, MST, and ST there is a reduction in YIELD compared to C (Control) and for T a significant increase can be seen.
Thurgem Dec 2023
thur.dec.23 |> summarise(m_y_t = mean(YIELD), .by = Code) |> ggplot(aes(x = Code, y = m_y_t)) + geom_bar(color = "black", alpha = 0.7, stat = "identity") + labs(title = "Bar chart of mean YIELD by Treaetment (Thurgem Dec 2023)", x = "Treatment", y = "Mean Yield") + theme_minimal()
For the dataset, it can be observed that for every treatment there is an increase in YIELD compared to treatment.
Kapsorok June 2022
kaps.jun.22 |> summarise(m_y_t = mean(YIELD), .by = Code) |> ggplot(aes(x = Code, y = m_y_t)) + geom_bar(color = "black", alpha = 0.7, stat = "identity") + labs(title = "Bar chart of mean YIELD by Treaetment (Kapsorok June 2022)", x = "Treatment", y = "Mean Yield") + theme_minimal()
For the dataset, it can be observed that for treatments M, MST and S there is a reduction in YIELD compared to C (Control) and for the rest of the treatments there is a significant increase in YIELD.
Visually tillage seem to have a negative effect on YIELD and manure combined with soil inoculum seems to have a slight positive effect on productivity but cannot be concluded.
Average MBC and MBN across treatment
Kabianga November 2023
It seem that manure slightly increase both MBC and MBN which is expected because MBC and MBN had a high corelation. It is also observed that tillage has a slightly negative effect.
kab.nov.22 |> summarise(m_mbc_t = mean(MBC), .by = Code) |> ggplot(aes(x = Code, y = m_mbc_t)) + geom_bar(color = "black", alpha = 0.7, stat = "identity") + labs(title = "Bar chart of mean MBC by Treaetment (Kabianga November 2022)", x = "Treatment", y = "Mean MBC") + theme_minimal()kab.nov.22 |> summarise(m_mbn_t = mean(MBN), .by = Code) |> ggplot(aes(x = Code, y = m_mbn_t)) + geom_bar(color = "black", alpha = 0.7, stat = "identity") + labs(title = "Bar chart of mean MBN by Treaetment (Kabianga November 2022)", x = "Treatment", y = "Mean MBN") + theme_minimal()
Primary analysis
Common function for brute force algorithm based on AIC. This was not implemented but only tested except for the dataset Kabianga November 2022 with dependent as YIELD because it yielded a sligtly better model. Hybrid stepwise algorithm is used based on AIC keeping the interaction effects as fixed for rest of the models. And for all models p<0.10 is considered significant.
#Pass the base model as the formulae and scope as a vector of variable names under consideration
bestModel <- function(base_model, scope, dataset)
{
all_combs <- list()
#seq_along creates n sequence of elements from the inputed vector
for (i in seq_along(scope))
{
#combn creates all possible n combinations, simplify is set to false to obtain the result as list
combs <- combn(scope, i, simplify = FALSE)
all_combs <- c(all_combs, combs)
}
combine <- function(all_combs)
{
#deparse convert type formula to string
paste(deparse(base_model), paste(all_combs, collapse =
" + "), sep = " + ")
}
#applies the function to list all_combs and returns a list which is converted to vector
formulas <- unlist(lapply(all_combs, combine))
fitted_model <- list()
for (formula in formulas)
{
model <- lm(as.formula(formula), data = dataset)
#adjust below to consider other criteria apart from AIC and also ajust the following logic accordingly after iteration
fitted_model[[formula]] <- list(models = model, AIC = AIC(model))
}
#takes in list, applies a function and returns vector
aic_values <- sapply(fitted_model, function(x) x$AIC)
#returns the fitted model and corresponding AIC
return(fitted_model[[which.min(aic_values)]])
}Kabianga November 2022
YIELD as the dependent variable.
im <- lm(YIELD ~ T * S * M, data = kab.nov.22) model_kab_Nov22_YIELD <- step(im, scope = list( lower = ~ T * S * M, upper = ~ T * S * M + Tsoil + Tair + MOIS + NH4 + NO3 + P + pH + EC + MBC + MBN + Nconc + Cconc ), direction = "both")Start: AIC=78.57 YIELD ~ T * S * M Df Sum of Sq RSS AIC + MOIS 1 25.6598 200.42 76.710 + MBN 1 15.9946 210.09 78.217 <none> 226.08 78.565 + NH4 1 11.2688 214.81 78.929 + pH 1 10.2982 215.78 79.073 + MBC 1 9.2920 216.79 79.222 + Tsoil 1 7.6671 218.41 79.461 + Tair 1 7.2481 218.83 79.523 + EC 1 6.5224 219.56 79.628 + NO3 1 5.6836 220.40 79.751 + P 1 2.5974 223.49 80.195 + Cconc 1 1.7285 224.35 80.320 + Nconc 1 0.0001 226.08 80.565 Step: AIC=76.71 YIELD ~ T + S + M + MOIS + T:S + T:M + S:M + T:S:M Df Sum of Sq RSS AIC + pH 1 13.4101 187.01 76.494 + MBN 1 13.0343 187.39 76.558 <none> 200.42 76.710 + P 1 10.0996 190.32 77.056 + Tair 1 9.3024 191.12 77.189 + NH4 1 9.1969 191.23 77.207 + MBC 1 6.5987 193.82 77.639 + Nconc 1 4.2286 196.19 78.028 + Tsoil 1 3.4636 196.96 78.152 - MOIS 1 25.6598 226.08 78.565 + Cconc 1 0.0500 200.37 78.702 + NO3 1 0.0310 200.39 78.705 + EC 1 0.0071 200.42 78.709 Step: AIC=76.49 YIELD ~ T + S + M + MOIS + pH + T:S + T:M + S:M + T:S:M Df Sum of Sq RSS AIC <none> 187.01 76.494 + Tair 1 11.2434 175.77 76.510 - pH 1 13.4101 200.42 76.710 + Tsoil 1 8.3602 178.65 77.031 + MBN 1 7.6768 179.34 77.153 + NH4 1 5.8166 181.20 77.483 + Nconc 1 4.1696 182.84 77.773 + P 1 4.1162 182.90 77.782 + MBC 1 2.3911 184.62 78.082 + Cconc 1 0.0688 186.94 78.482 + EC 1 0.0285 186.98 78.489 + NO3 1 0.0031 187.01 78.494 - MOIS 1 28.7718 215.78 79.073summary(model_kab_Nov22_YIELD)Call: lm(formula = YIELD ~ T + S + M + MOIS + pH + T:S + T:M + S:M + T:S:M, data = kab.nov.22) Residuals: Min 1Q Median 3Q Max -6.9018 -1.4359 0.3215 1.4152 7.5303 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 18.6138 14.5300 1.281 0.2135 T 1.2088 2.8672 0.422 0.6774 S 5.0615 2.2751 2.225 0.0367 * M 1.7076 2.0662 0.826 0.4174 MOIS 0.2926 0.1590 1.840 0.0793 . pH -3.4092 2.7143 -1.256 0.2223 T:S -3.2100 3.0478 -1.053 0.3037 T:M -2.1374 2.9200 -0.732 0.4719 S:M -6.6686 3.0928 -2.156 0.0423 * T:S:M 5.2884 4.3188 1.225 0.2337 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 2.916 on 22 degrees of freedom Multiple R-squared: 0.4955, Adjusted R-squared: 0.2891 F-statistic: 2.401 on 9 and 22 DF, p-value: 0.04534#implementation of brute force algorithm scope = c("Tsoil", "Tair", "MOIS", "NH4", "NO3", "P", "pH", "EC", "MBC", "MBN", "Nconc", "Cconc") base_model = YIELD ~ T * S * M #run only if necessar of running bcz it takes a quite a bit of time the model m below is the result of this function run. # m <- bestModel(base_model, scope, kab.nov.22) # summary(m$model) #model m below is the result of this function run (best_Model) m <- lm(YIELD ~ T * S * M + MOIS + NH4 + MBN + Nconc, data = kab.nov.22) summary(m)Call: lm(formula = YIELD ~ T * S * M + MOIS + NH4 + MBN + Nconc, data = kab.nov.22) Residuals: Min 1Q Median 3Q Max -4.3236 -1.3683 -0.1658 1.0606 4.9609 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -1.38995 8.65073 -0.161 0.87396 T 3.71723 2.81815 1.319 0.20206 S 6.29647 2.10707 2.988 0.00727 ** M 1.82889 1.95230 0.937 0.36004 MOIS 0.42629 0.16314 2.613 0.01665 * NH4 1.79159 1.49889 1.195 0.24596 MBN 0.12320 0.06815 1.808 0.08572 . Nconc -37.25304 15.70888 -2.371 0.02787 * T:S -2.37206 2.73333 -0.868 0.39578 T:M -2.22589 2.72349 -0.817 0.42339 S:M -7.41175 2.84954 -2.601 0.01709 * T:S:M 3.27673 3.82547 0.857 0.40184 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 2.695 on 20 degrees of freedom Multiple R-squared: 0.6081, Adjusted R-squared: 0.3926 F-statistic: 2.821 on 11 and 20 DF, p-value: 0.0212These results indicate that soil inoculum and the soil incolum combined with manure, moisture, MBN, and nitrogen concentration significantly influence yield, while other factors and interactions have less impact. Soil inoculum had a significant positive impact while soil incolum combined with manure had a negative impact.
MBC as the dependent variable
im <- lm(MBC ~ T * S * M, data = kab.nov.22) model_kab_Nov22_MBC <- step(im, scope = list( lower = ~ T * S * M, upper = ~ T * S * M + Tsoil + Tair + MOIS + NH4 + NO3 + P + pH + EC + MBN + Nconc + Cconc ), direction = "both")Start: AIC=261.84 MBC ~ T * S * M Df Sum of Sq RSS AIC + MBN 1 55169 14274 213.21 + Nconc 1 26477 42966 248.48 + Cconc 1 23970 45473 250.29 + NH4 1 16758 52685 255.00 + EC 1 14522 54921 256.33 + Tair 1 13595 55848 256.87 + Tsoil 1 13523 55920 256.91 + NO3 1 13448 55995 256.95 + pH 1 5741 63702 261.08 <none> 69443 261.84 + P 1 3725 65718 262.08 + MOIS 1 654 68789 263.54 Step: AIC=213.21 MBC ~ T + S + M + MBN + T:S + T:M + S:M + T:S:M Df Sum of Sq RSS AIC + P 1 2644 11629 208.66 + Cconc 1 1248 13026 212.29 <none> 14274 213.21 + NO3 1 819 13455 213.32 + EC 1 533 13740 214.00 + Nconc 1 519 13754 214.03 + pH 1 349 13925 214.42 + Tair 1 234 14040 214.69 + NH4 1 119 14154 214.94 + Tsoil 1 75 14199 215.05 + MOIS 1 49 14224 215.10 - MBN 1 55169 69443 261.84 Step: AIC=208.66 MBC ~ T + S + M + MBN + P + T:S + T:M + S:M + T:S:M Df Sum of Sq RSS AIC + pH 1 1704 9925 205.59 <none> 11629 208.66 + NH4 1 475 11154 209.32 + Cconc 1 312 11318 209.79 + Tair 1 68 11561 210.47 + MOIS 1 60 11569 210.49 + Nconc 1 22 11607 210.60 + Tsoil 1 7 11622 210.64 + EC 1 4 11625 210.65 + NO3 1 0 11629 210.66 - P 1 2644 14274 213.21 - MBN 1 54089 65718 262.08 Step: AIC=205.59 MBC ~ T + S + M + MBN + P + pH + T:S + T:M + S:M + T:S:M Df Sum of Sq RSS AIC + NH4 1 741 9184 205.10 <none> 9925 205.59 + Cconc 1 382 9543 206.33 + Tsoil 1 231 9694 206.83 + Tair 1 219 9706 206.87 + MOIS 1 58 9867 207.40 + Nconc 1 5 9920 207.57 + EC 1 3 9922 207.58 + NO3 1 1 9924 207.58 - pH 1 1704 11629 208.66 - P 1 4000 13925 214.42 - MBN 1 44873 54798 258.26 Step: AIC=205.1 MBC ~ T + S + M + MBN + P + pH + NH4 + T:S + T:M + S:M + T:S:M Df Sum of Sq RSS AIC <none> 9184 205.10 + Cconc 1 487 8697 205.36 - NH4 1 741 9925 205.59 + Tair 1 225 8959 206.31 + Tsoil 1 160 9024 206.54 + MOIS 1 74 9111 206.85 + Nconc 1 26 9159 207.01 + NO3 1 14 9170 207.06 + EC 1 10 9174 207.07 - pH 1 1970 11154 209.32 - P 1 4606 13790 216.11 - MBN 1 36389 45573 254.36summary(model_kab_Nov22_MBC)Call: lm(formula = MBC ~ T + S + M + MBN + P + pH + NH4 + T:S + T:M + S:M + T:S:M, data = kab.nov.22) Residuals: Min 1Q Median 3Q Max -29.70 -10.73 -3.20 10.10 40.76 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 368.5011 126.1290 2.922 0.00844 ** T -21.6695 16.0188 -1.353 0.19123 S 10.4362 16.8553 0.619 0.54280 M -4.2831 15.4943 -0.276 0.78505 MBN 4.3856 0.4927 8.902 2.15e-08 *** P -22.8369 7.2111 -3.167 0.00485 ** pH -46.2661 22.3393 -2.071 0.05151 . NH4 -14.6657 11.5453 -1.270 0.21856 T:S -36.6072 24.2079 -1.512 0.14612 T:M 2.4122 21.8744 0.110 0.91329 S:M -17.4707 23.1437 -0.755 0.45912 T:S:M 41.7015 33.6025 1.241 0.22896 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 21.43 on 20 degrees of freedom Multiple R-squared: 0.9063, Adjusted R-squared: 0.8548 F-statistic: 17.59 on 11 and 20 DF, p-value: 6.145e-08MBN as the dependent variable
im <- lm(MBN ~ T * S * M, data = kab.nov.22) model_kab_Nov22_MBN <- step(im, scope = list( lower = ~ T * S * M, upper = ~ T * S * M + Block + Tsoil + Tair + MOIS + NH4 + NO3 + P + pH + EC + MBC + Nconc + Cconc ), direction = "both")Start: AIC=161.23 MBN ~ T * S * M Df Sum of Sq RSS AIC + MBC 1 2378.00 615.24 112.60 + Nconc 1 1137.18 1856.06 147.94 + Block 1 1051.14 1942.10 149.38 + NH4 1 1037.41 1955.83 149.61 + Cconc 1 846.47 2146.77 152.59 + Tsoil 1 640.22 2353.02 155.53 + Tair 1 573.91 2419.34 156.42 + EC 1 538.12 2455.12 156.89 + NO3 1 434.67 2558.57 158.21 <none> 2993.24 161.23 + pH 1 180.42 2812.82 161.24 + MOIS 1 18.69 2974.55 163.03 + P 1 5.05 2988.19 163.17 Step: AIC=112.6 MBN ~ T + S + M + MBC + T:S + T:M + S:M + T:S:M Df Sum of Sq RSS AIC + NH4 1 89.79 525.45 109.55 + P 1 86.46 528.79 109.75 <none> 615.24 112.60 + Block 1 28.78 586.47 113.07 + Nconc 1 21.08 594.17 113.49 + Tsoil 1 17.77 597.47 113.66 + Tair 1 7.04 608.20 114.23 + EC 1 1.02 614.23 114.55 + NO3 1 0.46 614.78 114.58 + pH 1 0.38 614.87 114.58 + Cconc 1 0.30 614.94 114.58 + MOIS 1 0.17 615.08 114.59 - MBC 1 2378.00 2993.24 161.23 Step: AIC=109.55 MBN ~ T + S + M + MBC + NH4 + T:S + T:M + S:M + T:S:M Df Sum of Sq RSS AIC + P 1 103.39 422.07 104.54 <none> 525.45 109.55 + NO3 1 10.27 515.19 110.92 + Block 1 8.73 516.72 111.02 + Tsoil 1 4.64 520.81 111.27 + EC 1 2.46 522.99 111.40 + Tair 1 1.68 523.78 111.45 + Cconc 1 1.40 524.05 111.47 + Nconc 1 1.13 524.32 111.48 + pH 1 0.97 524.49 111.49 + MOIS 1 0.36 525.09 111.53 - NH4 1 89.79 615.24 112.60 - MBC 1 1430.38 1955.83 149.61 Step: AIC=104.54 MBN ~ T + S + M + MBC + NH4 + P + T:S + T:M + S:M + T:S:M Df Sum of Sq RSS AIC + pH 1 40.79 381.28 103.29 + Nconc 1 28.27 393.80 104.32 + Tair 1 26.83 395.24 104.44 <none> 422.07 104.54 + Block 1 19.79 402.27 105.00 + Tsoil 1 10.84 411.23 105.71 + EC 1 7.14 414.92 106.00 + MOIS 1 4.89 417.17 106.17 + NO3 1 3.66 418.40 106.26 + Cconc 1 1.65 420.41 106.42 - P 1 103.39 525.45 109.55 - NH4 1 106.72 528.79 109.75 - MBC 1 1518.73 1940.80 151.36 Step: AIC=103.29 MBN ~ T + S + M + MBC + NH4 + P + pH + T:S + T:M + S:M + T:S:M Df Sum of Sq RSS AIC + Tair 1 35.09 346.19 102.20 + Tsoil 1 24.36 356.91 103.18 + Block 1 23.93 357.35 103.22 <none> 381.28 103.29 + Nconc 1 17.90 363.38 103.75 - pH 1 40.79 422.07 104.54 + MOIS 1 4.98 376.30 104.87 + EC 1 4.71 376.57 104.89 + NO3 1 2.85 378.43 105.05 + Cconc 1 0.06 381.21 105.28 - NH4 1 117.87 499.14 109.91 - P 1 143.21 524.49 111.49 - MBC 1 1510.65 1891.93 152.55 Step: AIC=102.2 MBN ~ T + S + M + MBC + NH4 + P + pH + Tair + T:S + T:M + S:M + T:S:M Df Sum of Sq RSS AIC <none> 346.19 102.20 - Tair 1 35.09 381.28 103.29 + MOIS 1 7.22 338.97 103.53 + Nconc 1 3.01 343.18 103.92 + Tsoil 1 2.08 344.11 104.01 + Block 1 1.84 344.35 104.03 + NO3 1 1.74 344.45 104.04 + EC 1 1.38 344.81 104.07 + Cconc 1 0.21 345.98 104.18 - pH 1 49.05 395.24 104.44 - NH4 1 102.43 448.62 108.49 - P 1 176.78 522.97 113.40 - MBC 1 1297.49 1643.68 150.05summary(model_kab_Nov22_MBN)Call: lm(formula = MBN ~ T + S + M + MBC + NH4 + P + pH + Tair + T:S + T:M + S:M + T:S:M, data = kab.nov.22) Residuals: Min 1Q Median 3Q Max -6.6871 -2.5162 0.5216 2.5615 5.6404 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -39.0872 28.6689 -1.363 0.1887 T 2.3904 3.3107 0.722 0.4791 S -2.6106 3.4027 -0.767 0.4524 M 2.5024 3.0805 0.812 0.4267 MBC 0.1746 0.0207 8.439 7.52e-08 *** NH4 4.9858 2.1028 2.371 0.0285 * P 5.0085 1.6079 3.115 0.0057 ** pH 7.6942 4.6897 1.641 0.1173 Tair -0.2306 0.1661 -1.388 0.1813 T:S 8.2018 4.9141 1.669 0.1115 T:M -1.4562 4.3471 -0.335 0.7413 S:M 2.1094 4.6662 0.452 0.6563 T:S:M -7.2801 6.8728 -1.059 0.3028 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 4.269 on 19 degrees of freedom Multiple R-squared: 0.909, Adjusted R-squared: 0.8515 F-statistic: 15.81 on 12 and 19 DF, p-value: 2.073e-07MBC and MBN did not show any significant effects on treatment, but both had a very high significant effect on Phosphorus . MBC had a negative effect, whereas MBN had a positive effect on P.
Kabianga December 2023
im <- lm(YIELD ~ T * S * M, data = kab.dec.23) model_kab_Dec23_YIELD <- step(im, scope = list( lower = ~ T * S * M, upper = ~ T * S * M + P + BD ), direction = "both")Start: AIC=59.16 YIELD ~ T * S * M Df Sum of Sq RSS AIC <none> 123.28 59.159 + P 1 2.30009 120.98 60.556 + BD 1 0.36336 122.92 61.064summary(model_kab_Dec23_YIELD)Call: lm(formula = YIELD ~ T * S * M, data = kab.dec.23) Residuals: Min 1Q Median 3Q Max -3.6739 -1.2265 -0.0384 1.1020 3.6000 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 8.1330 1.1332 7.177 2.04e-07 *** T 0.6874 1.6026 0.429 0.6718 S 3.1477 1.6026 1.964 0.0612 . M -0.7659 1.6026 -0.478 0.6370 T:S -1.2792 2.2664 -0.564 0.5777 T:M 1.2228 2.2664 0.540 0.5945 S:M -1.2809 2.2664 -0.565 0.5772 T:S:M -2.4089 3.2052 -0.752 0.4596 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 2.266 on 24 degrees of freedom Multiple R-squared: 0.3125, Adjusted R-squared: 0.112 F-statistic: 1.559 on 7 and 24 DF, p-value: 0.1958Soil inoculum showed a positive effect, which was significant and soil inoculum combined with manure had a negative effect, even though not statistically significant. The effect of treatment seems comparable to the initial measurement, even though a reduction in YIELD is observed.
Thurgem June 2022
im <- lm(YIELD ~ T * S * M, data = thur.jun.22) model_Thur_Jun22 <- step(im, scope = list( lower = ~ T * S * M, upper = ~ T * S * M + Tair + Tsoil + MOIS + EC ), direction = "both")Start: AIC=26.79 YIELD ~ T * S * M Df Sum of Sq RSS AIC + MOIS 1 11.4629 33.370 19.341 <none> 44.832 26.790 + Tair 1 0.2215 44.611 28.632 + EC 1 0.1946 44.638 28.651 + Tsoil 1 0.1631 44.669 28.674 Step: AIC=19.34 YIELD ~ T + S + M + MOIS + T:S + T:M + S:M + T:S:M Df Sum of Sq RSS AIC <none> 33.370 19.341 + EC 1 0.4063 32.963 20.949 + Tair 1 0.1978 33.172 21.151 + Tsoil 1 0.0801 33.289 21.264 - MOIS 1 11.4629 44.832 26.790summary(model_Thur_Jun22)Call: lm(formula = YIELD ~ T + S + M + MOIS + T:S + T:M + S:M + T:S:M, data = thur.jun.22) Residuals: Min 1Q Median 3Q Max -2.4725 -0.3853 -0.0512 0.5106 2.5927 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 1.25689 1.23111 1.021 0.31790 T -0.70152 0.85285 -0.823 0.41920 S -0.28573 0.85199 -0.335 0.74039 M -1.33047 0.85187 -1.562 0.13198 MOIS 0.08097 0.02881 2.811 0.00992 ** T:S 0.90545 1.21536 0.745 0.46382 T:M 3.16417 1.20452 2.627 0.01507 * S:M 1.60516 1.20622 1.331 0.19632 T:S:M -3.33147 1.70502 -1.954 0.06298 . --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 1.205 on 23 degrees of freedom Multiple R-squared: 0.4221, Adjusted R-squared: 0.2211 F-statistic: 2.1 on 8 and 23 DF, p-value: 0.07843Tillage, when combined with manure, showed a significant positive effect. Moisture seems to have a very high significance with a p-value of 0.00992 with a positive effect.
Thurgem July 2023
Variable DW was avoided because of perfect correlation with YIELD.
im <- lm(YIELD ~ T * S * M, data = thur.july.23) model_Thur_Jul23 <- step(im, scope = list( lower = ~ T * S * M, upper = ~ T * S * M + Tair + Tsoil + MOIS + Nconc + Cconc + Block ), direction = "both")Start: AIC=120.51 YIELD ~ T * S * M Df Sum of Sq RSS AIC + Cconc 1 375.54 463.02 103.50 + Tair 1 107.83 730.74 118.11 <none> 838.56 120.51 + MOIS 1 14.73 823.83 121.94 + Tsoil 1 11.83 826.73 122.06 + Block 1 1.58 836.98 122.45 + Nconc 1 0.16 838.40 122.50 Step: AIC=103.51 YIELD ~ T + S + M + Cconc + T:S + T:M + S:M + T:S:M Df Sum of Sq RSS AIC <none> 463.02 103.50 + Tair 1 25.27 437.75 103.71 + Nconc 1 23.34 439.68 103.85 + MOIS 1 13.93 449.09 104.53 + Block 1 8.86 454.16 104.89 + Tsoil 1 0.89 462.13 105.44 - Cconc 1 375.54 838.56 120.51summary(model_Thur_Jul23)Call: lm(formula = YIELD ~ T + S + M + Cconc + T:S + T:M + S:M + T:S:M, data = thur.july.23) Residuals: Min 1Q Median 3Q Max -6.3913 -3.3870 0.1161 3.1192 7.5971 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -17.675 6.975 -2.534 0.018544 * T 9.336 3.417 2.732 0.011889 * S 3.218 3.236 0.995 0.330326 M 6.331 3.577 1.770 0.089949 . Cconc 21.942 5.080 4.319 0.000254 *** T:S -9.041 4.551 -1.987 0.059001 . T:M -9.553 4.874 -1.960 0.062238 . S:M -5.955 4.826 -1.234 0.229726 T:S:M 5.002 6.663 0.751 0.460427 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 4.487 on 23 degrees of freedom Multiple R-squared: 0.5609, Adjusted R-squared: 0.4082 F-statistic: 3.673 on 8 and 23 DF, p-value: 0.006751The treatment tillage had a significant positive effect, whereas tillage combined with soil inoculum and manure had a significant negative effect. Manure also showed a significant positive effect, although its effect is not as pronounced as tillage, with an estimate of 6.331, whereas tillage had 9.336. Carbon concentration had a very significance with a p-value of 0.000254 positive effect.
Thurgem Dec 2023
im <- lm(YIELD ~ T * S * M , data = thur.dec.23) model_Thur_dec23 <- step(im, scope = list( lower = ~ T * S * M, upper = ~ T * S * M + MOIS + P ), direction = "both")Start: AIC=10.88 YIELD ~ T * S * M Df Sum of Sq RSS AIC <none> 27.267 10.878 + MOIS 1 0.64431 26.622 12.113 + P 1 0.46632 26.800 12.326summary(model_Thur_dec23)Call: lm(formula = YIELD ~ T * S * M, data = thur.dec.23) Residuals: Min 1Q Median 3Q Max -1.6688 -0.6007 -0.1421 0.6171 2.6133 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 2.4113 0.5329 4.524 0.000139 *** T 1.0717 0.7537 1.422 0.167918 S 0.6397 0.7537 0.849 0.404415 M 0.8078 0.7537 1.072 0.294468 T:S -0.5609 1.0659 -0.526 0.603575 T:M -0.7722 1.0659 -0.724 0.475789 S:M 0.2662 1.0659 0.250 0.804918 T:S:M -0.5409 1.5074 -0.359 0.722870 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 1.066 on 24 degrees of freedom Multiple R-squared: 0.1976, Adjusted R-squared: -0.03645 F-statistic: 0.8443 on 7 and 24 DF, p-value: 0.5624Measurement done after a month (July 2023) showed no significant effects for any of the treatments, and the model obtained was fairly poor, with an R-squared of 19.76%.
Kapsorok June 2022
im <- lm(YIELD ~ T * S * M , data = kaps.jun.22) model_Kaps_jun22 <- step(im, scope = list( lower = ~ T * S * M, upper = ~ T * S * M + MOIS + Tsoil + Tair + EC + Block ), direction = "both")Start: AIC=51.3 YIELD ~ T * S * M Df Sum of Sq RSS AIC + Block 1 20.3776 76.072 45.710 <none> 96.449 51.305 + EC 1 5.3695 91.080 51.472 + MOIS 1 0.8644 95.585 53.017 + Tair 1 0.5081 95.941 53.136 + Tsoil 1 0.0257 96.423 53.296 Step: AIC=45.71 YIELD ~ T + S + M + Block + T:S + T:M + S:M + T:S:M Df Sum of Sq RSS AIC + MOIS 1 13.2794 62.792 41.571 <none> 76.072 45.710 + Tair 1 0.1182 75.953 47.660 + Tsoil 1 0.0763 75.995 47.678 + EC 1 0.0501 76.022 47.689 - Block 1 20.3776 96.449 51.305 Step: AIC=41.57 YIELD ~ T + S + M + Block + MOIS + T:S + T:M + S:M + T:S:M Df Sum of Sq RSS AIC <none> 62.792 41.571 + EC 1 0.585 62.207 43.271 + Tair 1 0.336 62.456 43.399 + Tsoil 1 0.242 62.551 43.448 - MOIS 1 13.279 76.072 45.710 - Block 1 32.793 95.585 53.017summary(model_Kaps_jun22)Call: lm(formula = YIELD ~ T + S + M + Block + MOIS + T:S + T:M + S:M + T:S:M, data = kaps.jun.22) Residuals: Min 1Q Median 3Q Max -2.9130 -0.5417 -0.0770 0.4879 4.2881 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 2.88600 1.34779 2.141 0.04358 * T 2.15809 1.19620 1.804 0.08492 . S -1.10692 1.19598 -0.926 0.36473 M -0.95527 1.19678 -0.798 0.43329 Block -1.35670 0.40025 -3.390 0.00264 ** MOIS 0.12675 0.05876 2.157 0.04220 * T:S 1.18085 1.69415 0.697 0.49309 T:M 0.06027 1.73371 0.035 0.97258 S:M 3.45978 1.71023 2.023 0.05539 . T:S:M -4.98098 2.49194 -1.999 0.05813 . --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 1.689 on 22 degrees of freedom Multiple R-squared: 0.594, Adjusted R-squared: 0.4279 F-statistic: 3.577 on 9 and 22 DF, p-value: 0.007063Measurements showed a significant negative effect for tillage combined with soil inoculum and manure, tillage showed a positive effect, and soi inoculum combined with manure showed a positive effect. Moisture showed a positive effect with a p-value of 0.04220.
Testing the significance of Block on dependent variable
summary(aov(YIELD ~ Block, data = kab.nov.22))Df Sum Sq Mean Sq F value Pr(>F) Block 1 9.0 8.978 0.745 0.395 Residuals 30 361.7 12.056summary(aov(YIELD ~ Block, data = kab.dec.23))Df Sum Sq Mean Sq F value Pr(>F) Block 1 1.8 1.801 0.304 0.585 Residuals 30 177.5 5.917summary(aov(YIELD ~ Block, data = thur.jun.22))Df Sum Sq Mean Sq F value Pr(>F) Block 1 2.58 2.576 1.401 0.246 Residuals 30 55.17 1.839summary(aov(YIELD ~ Block, data = thur.july.23))Df Sum Sq Mean Sq F value Pr(>F) Block 1 1.6 1.58 0.045 0.833 Residuals 30 1053.0 35.10summary(aov(YIELD ~ Block, data = thur.dec.23))Df Sum Sq Mean Sq F value Pr(>F) Block 1 1.43 1.426 1.314 0.261 Residuals 30 32.56 1.085summary(aov(YIELD ~ Block, data = kaps.jun.22))Df Sum Sq Mean Sq F value Pr(>F) Block 1 20.38 20.378 4.552 0.0412 * Residuals 30 134.29 4.476 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Except for Kapsorok Block did not have a significant effect on YIELD. The effect of Block for Kapsorok is significant with a p-value of 0.0412.
Creating table to summarize results for report.stargazer(model_Kaps_jun22, type = "text", title = "Comparison of Regression Models across datssets with YIELD as dependent variable")
Secondary analysis
Considering the variable count the datasets on Kabianga November 2022 was only chosen for principal component analysis.
Kabianga November 2023
Performing PCA and selecting principal components
pca_result.kab <- kab.nov.22 |>
select(-`Plot No.`, -M, -S, -T, - Treatment, -Code, -Block) |>
PCA(scale.unit = TRUE)fviz_eig(pca_result.kab, addlabels = TRUE, ylim = c(0, 50))Corelation table and biplot
var_coord <- pca_result.kab$var$coord[, 1:2]
selected_vars <- var_coord[c("Tsoil", "Tair", "MOIS", "YIELD", "NH4", "NO3", "P", "pH", "EC", "MBC", "MBN", "Nconc", "Cconc"), ]
correlation_table <- data.frame(Variable = rownames(selected_vars), round(selected_vars, 2))
corelation_tibble <- as_tibble(correlation_table) |>
rename(PC1 = Dim.1,
PC2 = Dim.2)
print(corelation_tibble)# A tibble: 13 × 3
Variable PC1 PC2
<chr> <dbl> <dbl>
1 Tsoil -0.61 0.52
2 Tair -0.69 0.5
3 MOIS 0.33 0.71
4 YIELD 0.2 0.74
5 NH4 0.67 0.18
6 NO3 0.72 -0.27
7 P -0.49 0.22
8 pH -0.02 -0.06
9 EC 0.76 -0.24
10 MBC 0.75 0.41
11 MBN 0.75 0.31
12 Nconc 0.88 -0.08
13 Cconc 0.8 0.19
K-means clustering based on principal components
Result may vary from report because set.seed(1234567) was added after report was written.
pc <- data.frame(pca_result.kab$ind$coord[, 1:2])
fviz_nbclust(pc, FUNcluster = kmeans, method = "wss")set.seed(12334567)
k_means <- kmeans(pc, centers = 4)
fviz_pca_biplot(pca_result.kab,
habillage = as.factor(k_means$cluster),
addEllipses = TRUE
)fviz_pca_biplot(pca_result.kab,
habillage = as.factor(kab.nov.22$Code)
)p_cat <- tibble(
Code = kab.nov.22$Code,
Cluster = as.factor(k_means$cluster)
)
ggplot(p_cat, aes(x = Code, fill=Cluster)) +
geom_bar(stat = "count") +
theme_minimal()It seems as though manure treatment is associated with improved soil properties which are positively correlated with YIELD, except for EC and NO3, which seem to be poorly correlated. Tilage seems to have a very little to no impact on YIELD.
Just a plot to visualise variables and their correlation with PC1 and PC2
correlation_data <- data.frame(
Variable = rownames(var_coord),
PC1 = pca_result.kab$var$coord[, 1],
PC2 = pca_result.kab$var$coord[, 2]
)
# Correlation plot
ggplot(correlation_data, aes(x = PC1, y = PC2, label = Variable)) +
geom_point(size = 3, color = "blue") +
geom_text(vjust = 1.5, color = "black") +
theme_minimal() +
labs(title = "Correlations of Variables with Principal Components",
x = "PC1 (Principal Component 1)",
y = "PC2 (Principal Component 2)")Assesing the effect of solar radiation flux and precipitation
Creating tibbble for analysis
wi_kab_nov_22 <- read_csv("data/Weather data/wi_kab_dec_23.csv")
wi_kab_dec_23 <- read_csv("data/Weather data/wi_kab_nov_22.csv")
wi_thur_jun_22 <- read_csv("data/Weather data/wi_thur_jun_22.csv")
wi_thur_jul_23 <- read_csv("data/Weather data/wi_thur_jul_23.csv")
wi_thur_dec_23 <- read_csv("data/Weather data/wi_thur_dec_23.csv")
wi_kaps_jun_22 <- read_csv("data/Weather data/wi_kaps_jun_22.csv")
a <- kab.nov.22 |>
select(YIELD, Block, Code, MOIS) |>
summarise(mean_yield = mean(YIELD), .by = c(Code)) |>
mutate(Mean_precip = mean(wi_kab_nov_22$Precipitation),
Mean_srf = mean(wi_kab_nov_22$solar_radiation_flux),
date = "nov_22")
b <- kab.dec.23 |>
select(YIELD, Block, Code) |>
summarise(mean_yield = mean(YIELD), .by = c(Code)) |>
mutate(Mean_precip = mean(wi_kab_dec_23$Precipitation),
Mean_srf = mean(wi_kab_dec_23$solar_radiation_flux),
date = "dec_23")
c <- thur.jun.22 |>
select(YIELD, Block, Code, MOIS) |>
summarise(mean_yield = mean(YIELD), .by = c(Code)) |>
mutate(Mean_precip = mean(wi_thur_jun_22$Precipitation),
Mean_srf = mean(wi_thur_jun_22$solar_radiation_flux),
date = "jun_22")
d <- thur.july.23 |>
select(YIELD, Block, Code) |>
summarise(mean_yield = mean(YIELD), .by = c(Code)) |>
mutate(Mean_precip = mean(wi_thur_jul_23$Precipitation),
Mean_srf = mean(wi_thur_jul_23$solar_radiation_flux),
date = "jul_23")
e <- thur.dec.23 |>
select(YIELD, Block, Code) |>
summarise(mean_yield = mean(YIELD), .by = c(Code)) |>
mutate(Mean_precip = mean(wi_thur_dec_23$Precipitation),
Mean_srf = mean(wi_thur_dec_23$solar_radiation_flux),
date = "dec_23")
weather_info_kab <- bind_rows(a, b)
weather_info_thur <- bind_rows(c, d ,e)
weather_info_kaps <- kaps.jun.22 |>
select(YIELD, Block, Code) |>
summarise(mean_yield = mean(YIELD), .by = c(Code)) |>
mutate(Mean_precip = mean(wi_kaps_jun_22$Precipitation),
Mean_srf = mean(wi_kaps_jun_22$solar_radiation_flux),
date = "jun_22")
rm(a,b,c,d,e)Kabianga
ggplot(weather_info_kab, aes(x = Mean_precip, y = mean_yield, group = Code, color = Code)) + geom_line(aes(linetype = Code)) + geom_point() + geom_text_repel(aes(label = date)) + facet_wrap(~Code) + theme_minimal() + labs(title = "Yield Trends Over Time by Treatment For Precipitation", x = "Mean precipitation", y = "Yield")ggplot(weather_info_kab, aes(x = log(Mean_srf), y = mean_yield, group = Code, color = Code)) + geom_line(aes(linetype = Code)) + geom_point() + geom_text_repel(aes(label = date)) + facet_wrap(~Code) + theme_minimal() + labs(title = "Yield Trends Over Time by Treatment For Solar radiation Flux ", x = "Mean solar radiation flux (log scale)", y = "Yield")summary(lm(mean_yield ~ Mean_precip:Code, data = weather_info_kab))Call: lm(formula = mean_yield ~ Mean_precip:Code, data = weather_info_kab) Residuals: Min 1Q Median 3Q Max -1.59945 -1.16581 0.05666 0.87641 2.26297 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 35.7426 2.5642 13.939 2.31e-06 *** Mean_precip:CodeC -5.3420 0.6415 -8.328 7.05e-05 *** Mean_precip:CodeM -5.2672 0.6415 -8.211 7.71e-05 *** Mean_precip:CodeMS -5.1600 0.6415 -8.044 8.80e-05 *** Mean_precip:CodeMST -5.6693 0.6415 -8.838 4.80e-05 *** Mean_precip:CodeMT -5.4248 0.6415 -8.457 6.38e-05 *** Mean_precip:CodeS -4.5853 0.6415 -7.148 0.000186 *** Mean_precip:CodeST -5.1048 0.6415 -7.958 9.43e-05 *** Mean_precip:CodeT -5.4642 0.6415 -8.518 6.09e-05 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 1.769 on 7 degrees of freedom Multiple R-squared: 0.9279, Adjusted R-squared: 0.8454 F-statistic: 11.25 on 8 and 7 DF, p-value: 0.002275summary(lm(mean_yield ~ Mean_srf:Code, data = weather_info_kab))Call: lm(formula = mean_yield ~ Mean_srf:Code, data = weather_info_kab) Residuals: Min 1Q Median 3Q Max -2.0193 -0.7584 -0.0061 0.7237 1.8153 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -6.176e+01 7.236e+00 -8.535 6.02e-05 *** Mean_srf:CodeC 3.680e-06 3.604e-07 10.211 1.86e-05 *** Mean_srf:CodeM 3.711e-06 3.604e-07 10.298 1.76e-05 *** Mean_srf:CodeMS 3.715e-06 3.604e-07 10.311 1.75e-05 *** Mean_srf:CodeMST 3.597e-06 3.604e-07 9.982 2.17e-05 *** Mean_srf:CodeMT 3.641e-06 3.604e-07 10.103 2.00e-05 *** Mean_srf:CodeS 3.849e-06 3.604e-07 10.681 1.38e-05 *** Mean_srf:CodeST 3.711e-06 3.604e-07 10.299 1.76e-05 *** Mean_srf:CodeT 3.636e-06 3.604e-07 10.090 2.02e-05 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 1.538 on 7 degrees of freedom Multiple R-squared: 0.9454, Adjusted R-squared: 0.8831 F-statistic: 15.16 on 8 and 7 DF, p-value: 0.0008942Thurgum
ggplot(weather_info_thur, aes(x = Mean_precip, y = mean_yield, group = Code, color = Code)) + geom_line(aes(linetype = Code)) + geom_point() + geom_text_repel(aes(label = date)) + facet_wrap(~Code) + theme_minimal() + labs(title = "Yield Trends Over Time by Treatment For Precipitation", x = "Mean precipitation", y = "Yield")ggplot(weather_info_thur, aes(x = log(Mean_srf), y = mean_yield, group = Code, color = Code)) + geom_line(aes(linetype = Code)) + geom_point() + geom_text_repel(aes(label = date)) + facet_wrap(~Code) + theme_minimal() + labs(title = "Yield Trends Over Time by Treatment For Solar radiation Flux ", x = "Mean solar radiation flux (log scale)", y = "Yield")summary(lm(mean_yield ~ Mean_precip:Code, data = weather_info_thur))Call: lm(formula = mean_yield ~ Mean_precip:Code, data = weather_info_thur) Residuals: Min 1Q Median 3Q Max -4.2480 -2.9100 -0.0426 2.0472 6.8268 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 9.370 1.516 6.180 1.76e-05 *** Mean_precip:CodeC -2.885 1.716 -1.681 0.114 Mean_precip:CodeM -2.837 1.716 -1.653 0.119 Mean_precip:CodeMS -2.164 1.716 -1.261 0.227 Mean_precip:CodeMST -3.177 1.716 -1.851 0.084 . Mean_precip:CodeMT -2.216 1.716 -1.291 0.216 Mean_precip:CodeS -2.639 1.716 -1.537 0.145 Mean_precip:CodeST -2.689 1.716 -1.567 0.138 Mean_precip:CodeT -2.223 1.716 -1.295 0.215 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 3.887 on 15 degrees of freedom Multiple R-squared: 0.3227, Adjusted R-squared: -0.03853 F-statistic: 0.8933 on 8 and 15 DF, p-value: 0.5452summary(lm(mean_yield ~ Mean_srf:Code, data = weather_info_thur))Call: lm(formula = mean_yield ~ Mean_srf:Code, data = weather_info_thur) Residuals: Min 1Q Median 3Q Max -4.5003 -2.4308 0.3415 1.9385 4.8084 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 7.197e+01 1.566e+01 4.595 0.000350 *** Mean_srf:CodeC -3.230e-06 7.686e-07 -4.203 0.000769 *** Mean_srf:CodeM -3.251e-06 7.686e-07 -4.229 0.000728 *** Mean_srf:CodeMS -3.182e-06 7.686e-07 -4.140 0.000873 *** Mean_srf:CodeMST -3.310e-06 7.686e-07 -4.307 0.000623 *** Mean_srf:CodeMT -3.169e-06 7.686e-07 -4.123 0.000903 *** Mean_srf:CodeS -3.215e-06 7.686e-07 -4.183 0.000799 *** Mean_srf:CodeST -3.242e-06 7.686e-07 -4.218 0.000745 *** Mean_srf:CodeT -3.165e-06 7.686e-07 -4.118 0.000913 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 3.093 on 15 degrees of freedom Multiple R-squared: 0.571, Adjusted R-squared: 0.3422 F-statistic: 2.495 on 8 and 15 DF, p-value: 0.06039
In Kabianga, productivity was lower during the months with high rainfall. However, when considering solar radiation flux, yield increased with higher levels of solar radiation. In contrast, at Thurgem, YIELD decreased as both rainfall and solar radiation flux increased.